SUBROUTINE vfunapprox_theta_thetav_ability_age(agenin,agemayorin,emaxin,betain)
  USE Commonvars
  IMPLICIT NONE
  INCLUDE 'mpif.h'

  INTEGER, INTENT(IN)  :: agenin, agemayorin(agenin)
  REAL(8), INTENT(IN)  :: emaxin(Ntheta,Ntheta_vot,Nab,Nwealth,agenin)
  REAL(8), INTENT(OUT) :: betain(MMpmnew_age)
  REAL(8), ALLOCATABLE :: y(:), xvar(:,:), res(:,:,:,:,:), vfunout(:,:,:,:,:), ones(:), xpx(:,:), xpy(:), invxx(:,:), res_final(:), work(:), betastart(:)
  REAL(8)              :: Rsquared, dble_agemayorin(agenin), av_res(Nwealth), av_res_temp(Nwealth), v_slope, v_slope_p
  INTEGER              :: ii, i, j, iwealth, ia, itheta, ithetai, info, MM_t, agei, l_wealth, lwork
  CHARACTER :: uplo='U'

  ii = Ntheta*Ntheta_vot*Nab*Nwealth*agenin
  dble_agemayorin(:) = DBLE(agemayorin(:))

  ALLOCATE(y(ii), res_final(ii), res(Ntheta,Ntheta_vot,Nab,agenin,Nwealth), ones(ii), vfunout(Ntheta,Ntheta_vot,Nab,agenin,Nwealth))
  y = 0d0; res = 0d0

  MM_t = Npol_age+Npol_ab+Npol_theta+Npol_theta_vot
  ALLOCATE(xvar(ii,MM_t),xpx(MM_t,MM_t), xpy(MM_t), invxx(MM_t,MM_t), betastart(MM_t))

  xpx = 0d0; xpy = 0d0
  betastart = 0d0

  !Average the residuals for each iwealth
  av_res_temp = 0d0
  DO iwealth = 1, Nwealth
     DO agei = 1, agenin
        DO ia = 1, Nab
           DO itheta = 1, Ntheta_vot
              DO ithetai = 1, Ntheta
                 av_res_temp(iwealth) = av_res_temp(iwealth) + emaxin(ithetai,itheta,ia,iwealth,agei)
              END DO
           END DO
        END DO
     END DO
  END DO
  av_res_temp = av_res_temp/DBLE(Ntheta*Ntheta_vot*Nab*agenin)

  !Compute the residuals V - (alpha_1+alpha_2*q+alpha_3*q^2)
  DO iwealth = 1, Nwealth
     DO agei = 1, agenin
        DO ia = 1, Nab
           DO itheta = 1, Ntheta_vot
              DO ithetai = 1, Ntheta
                 res(ithetai,itheta,ia,agei,iwealth) = emaxin(ithetai,itheta,ia,iwealth,agei) - av_res_temp(iwealth)
              END DO
           END DO
        END DO
     END DO
  END DO

  !We compute the coefficients on public consumption in V = f(sav)+alpha_1+alpha_2*q+alpha_3*q^2
  ones = 1d0
  i = 0
  DO iwealth = 1, Nwealth
     DO agei = 1, agenin
        DO ia = 1, Nab
           DO itheta = 1, Ntheta_vot
              DO ithetai = 1, Ntheta
                 
                 i = i + 1
                 y(i) = res(ithetai,itheta,ia,agei,iwealth)
                 xvar(i,1:MM_t) = (/ ones(i), a_discr(ia), a_discr(ia)**2d0, dble_agemayorin(agei), dble_agemayorin(agei)**2d0, theta_vot(itheta), theta(ithetai) /)

              END DO
           END DO
        END DO
     END DO
  END DO
  
  xpx = MATMUL(TRANSPOSE(xvar(:,:)),xvar(:,:))/SIZE(y)
  xpy = MATMUL(TRANSPOSE(xvar(:,:)),y(:))/SIZE(y)
     
  invxx = 0d0
  DO j=1, MM_t
     invxx(j,j)=1d0
  END DO
  
  LWORK = MM_t
  IF (ALLOCATED(WORK)) DEALLOCATE(WORK)
  ALLOCATE(WORK(LWORK))
  CALL DPOSV(uplo,MM_t,MM_t,xpx,MM_t,invxx,MM_t,info)
  IF (info /= 0) THEN
     DO j = 1, ii
        WRITE(*,'(a,1i3,7f12.6)') 'y,x',j,y(j),xvar(j,1:MM_t)
     ENDDO
     DO iwealth = 1, Nwealth
        DO agei = 1, agenin
           DO ia = 1, Nab
              DO itheta = 1, Ntheta_vot
                 DO ithetai = 1, Ntheta
                    
                    WRITE(*,'(a,5i4,2f12.6)') 'emaxin', ithetai,itheta,ia,agei,iwealth, emaxin(ithetai,itheta,ia,iwealth,agei), av_res_temp(iwealth)

                 END DO
              END DO
           END DO
        END DO
     END DO
     WRITE(*,*) 'MATRIX COULD NOT INVERT vfun_age 1', 'info', info
  END IF
  
  betastart(:) = MATMUL(invxx,xpy)

  IF (flag_Tend == 1) THEN
     v_slope_p = (av_res_temp(2)-av_res_temp(1))/(wealth(2) - wealth(1))
     av_res(Nwealth) = av_res_temp(Nwealth)
     DO iwealth = Nwealth-1, 1, -1
        IF (av_res_temp(iwealth+1) - av_res_temp(iwealth) > small_no) THEN
           av_res(iwealth) = av_res_temp(iwealth)
           l_wealth = iwealth
           V_slope = (av_res_temp(iwealth+1)-av_res_temp(iwealth))/(wealth(iwealth+1) - wealth(iwealth))
        ELSE        
           V_slope = inc_slope*V_slope_p+V_slope_p
           av_res(iwealth) = V_slope*(wealth(iwealth) - wealth(iwealth+1)) + av_res(iwealth+1)

           IF (flag_here == 9823) THEN
              IF (av_res(iwealth+1) - av_res(iwealth) < 0d0) THEN
                 write(*,'(a,1i4,5f15.7)') 'beta vfunapprox_age 1', iwealth, av_res(iwealth), V_slope, V_slope_p, av_res(iwealth+1), wealth(iwealth) - wealth(iwealth+1)              
                 write(*,'(a,10f15.7)') 'beta vfunapprox_age 2', av_res(:)
              END IF
           END IF

        END IF
        V_slope_p = V_slope
     END DO

!!$     IF (flag_here == 9823) write(*,'(a,10f15.7)') 'beta vfunapprox_age', av_res(:)

  ELSE
     av_res = av_res_temp
  END IF
  betain(:) = (/ av_res(:), betastart(:) /)

  !Compare the actual value function with the approximation
  i = 0
  DO iwealth = 1, Nwealth
     DO agei = 1, agenin
        DO ia = 1, Nab
           DO itheta = 1, Ntheta_vot
              DO ithetai = 1, Ntheta
                 
                 i = i + 1
                 vfunout(ithetai,itheta,ia,agei,iwealth) = av_res(iwealth) + betastart(1) + betastart(2)*a_discr(ia) + betastart(3)*a_discr(ia)**2d0 + betastart(4)*dble_agemayorin(agei)+betastart(5)*dble_agemayorin(agei)**2d0 + betastart(6)*theta_vot(itheta) + betastart(7)*theta(ithetai)
                 res_final(i) = emaxin(ithetai,itheta,ia,iwealth,agei) - vfunout(ithetai,itheta,ia,agei,iwealth)
                 y(i) = emaxin(ithetai,itheta,ia,iwealth,agei)

!!$           IF (myrank == 0) write(*,'(a,4i4,6f20.10)') 'vfun_age', iq,agei,ia,iwealth,vfunout(ia,agei,iwealth),emaxin(ithetai,itheta,ia,iwealth,agei),&
!!$                av_res_temp(iwealth) + betastart(1) + betastart(2)*a_discr(ia) + betastart(3)*a_discr(ia)**2d0 + betastart(5)*qconsvec(iq)**2d0 + betastart(6)*dble_agemayorin(agei)+betastart(7)*dble_agemayorin(agei)**2d0, av_res(iwealth), betastart(1) + betastart(2)*a_discr(ia) + betastart(3)*a_discr(ia)**2d0 + betastart(6)*dble_agemayorin(agei)+betastart(7)*dble_agemayorin(agei)**2d0

              END DO
           END DO
        END DO
     END DO
  END DO

  Rsquared = (SUM((y-SUM(y)/DBLE(SIZE(y)))**2d0) - SUM(res_final**2d0)) / SUM((y-SUM(y)/DBLE(SIZE(y)))**2d0)

  IF (Rsquared < 0d0 .AND. flag_Tend == 0) THEN
!!$  IF (Rsquared < 0.75d0) THEN
!!$  IF (flag_here == 1) THEN
     write(*,*) 'Rsquared vfunapprox_age', Rsquared
     DO agei = 1, agenin           
        DO ia = 1, Nab
           DO iwealth = 1, Nwealth
              DO itheta = 1, Ntheta_vot
                 DO ithetai = 1, Ntheta

                    write(*,'(a,5i4,5f20.10)') 'vfunapprox_theta_thetav_ability_age', ithetai,itheta,agei,ia,iwealth,emaxin(ithetai,itheta,ia,iwealth,agei), vfunout(ithetai,itheta,ia,agei,iwealth), av_res(iwealth) + betastart(1) + betastart(2)*a_discr(ia) + betastart(3)*a_discr(ia)**2d0 + betastart(4)*dble_agemayorin(agei)+betastart(5)*dble_agemayorin(agei)**2d0 + betastart(6)*theta_vot(itheta) + betastart(7)*theta(ithetai), av_res_temp(iwealth), betastart(1) + betastart(2)*a_discr(ia) + betastart(3)*a_discr(ia)**2d0 + betastart(4)*dble_agemayorin(agei)+betastart(5)*dble_agemayorin(agei)**2d0 + betastart(6)*theta_vot(itheta) + betastart(7)*theta(ithetai)

                 END DO
              END DO           
           END DO
        END DO
     END DO
  END IF

!!$  IF (myrank == 0)write(*,*) 'betain', betain
!!$  IF (myrank == 0)write(*,'(a,10f20.10)') 'av_res', av_res
!!$  IF (myrank == 0)write(*,'(a,10f20.10)') 'av_res_temp', av_res_temp
!!$  IF (myrank == 0)write(*,*) ''


  DEALLOCATE(y,xvar,res,vfunout)

END SUBROUTINE vfunapprox_theta_thetav_ability_age

